#include <iostream>
#include <fstream>
#include <limits>
#include <vector>
#include <sstream>
#include <algorithm>
#include <boost/regex.hpp>
#include <math.h>
#include <assert.h>

using namespace std;
using namespace boost;

const double xmin = -2.0;
const double xmax = 50.0;
const double xstep = 0.01;
const double bandwidth = 0.5;
double triggerlevel = 1.0;
const int numevalevents = 10000;

vector<double> x;

int processHeader(istream *infile)
{
    int created = 0;
    
    static const regex re(".*producing ([0-9]+) optical photons.*");
    
    smatch mat;
    
    // skip comments and newlines, add created photon numbers
    while (*infile) {
        string commentline;
        switch (infile->peek()) {
        case '\n':
            infile->get();
            break;
        case '#':
            getline(*infile, commentline);
            if (commentline.compare(0, 11, "# Time (ns)") == 0) {
                return created;
            }
            else if (commentline.compare(0, 17, "# scintillating: ") == 0) {
                regex_match(commentline, mat, re);
                if (mat.size() == 2) {
                    std::string m(mat[1].str());
                    int c = 0;
                    stringstream(m) >> c;
                    created += c;
                }
            }
            break;
        default:
            return created;
        }
    }
    return created;
}

double readSingleValue(istream *infile)
{
    double val = 0.0;
    *infile >> val;
    infile->ignore(numeric_limits<int>::max(), '\n');
    return val;
}

inline double kernel(double x) {
    static const double a = 35.0/32.0;

    const double f = 1.0 - (x * x);
    const double result = a * f * f * f;
    assert(result >= 0.0);
    return result;
}

bool abovetrigger(const double &a)
{
    return a >= triggerlevel;
}


void processEvent(const vector<double> &event, double *maximum, double *trigger) {
    // determine scaling factor for canonical triweight kernel
    static const double h = pow(9450.0/143.0, 0.2) * bandwidth;
    // bins to update for a single photon
    static const double halfwindow = h;

    // initialize density result vector
    vector<double> y(x.size(), 0.0);

    // add bump for each event
    for (vector<double>::const_iterator e=event.begin(); e!=event.end(); e++) {
        // determine bin of event (lower index)
        int minbin = ceil(((*e-halfwindow) - xmin) / xstep);
        minbin = max(0, minbin);
        int maxbin = floor(((*e+halfwindow) - xmin) / xstep);
        maxbin = min(int(y.size())-1, maxbin);
        for (int i=minbin; i<=maxbin; i++) {
            double evalat = (x.at(i)-(*e))/h;
            y[i] += kernel(evalat) / h;
        }
    }

    // determine max and increment histogram bin
    *maximum = x.at(max_element(y.begin(), y.end())-y.begin());

    // determine trigger position and increment histogram bin
    vector<double>::const_iterator t = find_if(y.begin(), y.end(), &abovetrigger);
    if (t != y.end())
        *trigger = x.at(t-y.begin());
    else
        *trigger = numeric_limits<double>::quiet_NaN();
}


int main(int argc, char* argv[])
{

    if (argc < 2) {
        cout << "Need at least one argument! Usage: kerneldens basename" << endl;
        return 1;
    }

    string outfilename(argv[1]);
    outfilename.append("-density.dat");
    {
        ifstream testoutfile(outfilename.c_str());
        if (testoutfile.good()) {
            cout << "Outfile already exists. Please delete/rename before running kerneldens again!" << endl;
            return 3;
        }
        testoutfile.close();
    }

    cout << "Using basename " << argv[1] << endl;

    // setting triggerlevel (max signal of 10 coincident photons)
    triggerlevel = 10.0 * (35.0/32.0 / (pow(9450.0/143.0, 0.2) * bandwidth));

    // initializing x
    for (int i=0; i<=int(ceil((xmax-xmin)/xstep)+0.5); i++)
        x.push_back(xmin + double(i)* xstep);

    // initializing maximum histogram
    //vector<int> maxhist(x.size(), 0);
    // initializing trigger histogram
    //vector<int> triggerhist(x.size(), 0);
    vector<int> created(numevalevents, 0);
    vector<double> maximum(numevalevents, numeric_limits<double>::quiet_NaN());
    vector<double> trigger(numevalevents, numeric_limits<double>::quiet_NaN());

    int numevents = 0;
    int fileno = 0;
    int goodevents = 0;

    while (numevents < numevalevents) {

        ostringstream filename;
        filename << argv[1] << "." << ++fileno << ".dat";
        ifstream infile;
        infile.open(filename.str().c_str());

        if (!infile) {
            cout << "Cannot open input file no " << fileno << "!" << endl;
            cout << "Filename was " << filename.str() << endl;
            break;
        }
        else
            cout << "Opened file no " << fileno << "." << endl;

        while (numevents < numevalevents) { // is also stopped on eof

            int c = processHeader(&infile);

            if (infile.eof())
                break;

            // read current event
            vector<double> event;
            double val = readSingleValue(&infile);
            while (infile) {
                event.push_back(val);
                val = readSingleValue(&infile);
            }

            // recover from failed double-read
            infile.clear();

            double m = numeric_limits<double>::quiet_NaN(), t = numeric_limits<double>::quiet_NaN();
            if (!event.empty() && c>0) {
                processEvent(event, &m, &t);
                created[goodevents] = c;
                maximum[goodevents] = m;
                trigger[goodevents] = t;
                goodevents++;
            }

            numevents++;

        }

	infile.close();

    }

    cout << "Evaluated a total of " << numevents << " events." << endl;

    if (numevents == 0)
    	return 2;

    ofstream outfile;
    outfile.open(outfilename.c_str());

    if (!outfile.is_open()) {
        cout << "Cannot open output file!" << endl;
        return 2;
    }

    outfile << "# Kernel density evaluation of " << argv[1] << ".?.dat" << endl;
    outfile << "# Evaluated a total of " << numevents << " events." << endl;
    outfile << "# Chosen bandwidth: " << bandwidth << endl;
    outfile << "# Chosen trigger level: " << triggerlevel << endl << endl;
    outfile << "\"Created Photons\"\t\"Time of Maximum\"\t\"Time of Trigger\"" << endl;
    for (int i=0; i<goodevents; i++)
        outfile << created.at(i) << "\t" << maximum.at(i) << "\t" << trigger.at(i) << endl;

    if (numevents < numevalevents)
        return 1;

    return 0;
}

